function r = gauss_int1d(function_f, a, b)
gauss_points_1d = [-0.9681602395076260898355762, -0.8360311073266357942994298, -0.6133714327005903973087020, -0.3242534234038089290385380, 0, 0.3242534234038089290385380, 0.6133714327005903973087020, 0.8360311073266357942994298, 0.9681602395076260898355762];
gauss_weights_1d = [0.0812743883615744119718922, 0.1806481606948574040584720, 0.2606106964029354623187429, 0.3123470770400028400686304, 0.3302393550012597631645251, 0.3123470770400028400686304, 0.2606106964029354623187429, 0.1806481606948574040584720, 0.0812743883615744119718922];
A = 0.5*(b-a);
B = 0.5*(a+b);
local_points = A*gauss_points_1d + B;
local_weights = A*gauss_weights_1d;
r = 0;
for i = 1:9
    r = r + local_weights(i) * function_f(local_points(i));
end
end